---
title: "Integrated ecosystem assessment on MFA (multiple factor analysis) for Black sea pelagic ecosystem, 1994-2017"
author: "Author: Piatinskii M., Azov-black sea branch of VNIRO"
date: 'Report build date: `r Sys.time()`'
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 1. About
MFA - multiple Factor Analysis that studies several groups of variables defined on the same set of individuals. The core of the method is a factor analysis applied to the whole set of variables in which each group of variables is weighted. This point of view leads to a representation of individuals and variables, as in any factor analysis. Owing to the weighting, this factor analysis can be interpreted as a canonical analysis. This point of view induces a display in which representations of the set of individuals associated to each group of variables are superposed (these displays are akin to that of procrustes analysis).

MFA apply serveral methods based on data type:
  - for the quantitative sets, to the one used in PCA; 
  - for the categorical sets, to the one used in MCA; 
  - for the frequency sets, to the one used in MFACT.

MFA can deal with NA in data by avaraging standartization procedure. Quantitative data can be standartized by robust procedure.

Refs:

  - Escofier, B. and Pages, J. (1994) Multiple Factor Analysis (AFMULT package). Computational Statistics and Data Analysis, 18, 121-140.
Becue-Bertaut, M. and Pages, J. (2008) 
  - Multiple factor analysis and clustering of a mixture of quantitative, categorical and frequency data. Computational Statistice and Data Analysis, 52, 3255-3268.
  

## 2. Input data
Here you can find input data to PCA analysis.

```{r echo=FALSE}
knitr::kable(data, digits = 3, caption = "Table 2.1. PCA input data")
```

## 3. Data diagnostics
There is preliminary data diagnostics procedure shown. Anomalies and traffic light plot should be reviewed to make primary data identification. 

### 3.1. Anomalies
This section shows the anomaly plots for all input data columns. Anomalies (a{x}) measurement calculated in simple way:

$$
\overline{x}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}={\frac {x_{1}+x_{2}+\cdots +x_{n}}{n}} \\
a_{x}=x_{i}-\overline{x}
$$

where x - observations variants row. 

```{r echo=FALSE, fig.cap="Fig. 3.1.1. Environment anomaly diagnostcs"}
knitr::include_graphics("output/anomaly/env.png")
```

```{r echo=FALSE, fig.cap="Fig. 3.1.2. Sprat and plankton anomaly diagnostics"}
knitr::include_graphics("output/anomaly/plankton-sprat.png")
```

### 3.2. Traffic light
Traffic light plot show input data continious changes based on quantile deviatians. Quantile calculated for sorted data in next intervals: 

  - 0 ... 0.2 - shows high significant reduction in observation year by factor
  - 0.2 ... 0.4 - shows moderate reduction in observation year by factor
  - 0.4 ... 0.6 - shows no significant changes in observation year by factor
  - 0.6 ... 0.8 - shows moderate increase in observation year by factor
  - 0.8 ... 1 - show high significant increase in observation year by factor
  
```{r echo=FALSE, fig.cap = "3.2.1. Traffic light diagnostics plot"}
knitr::include_graphics("output/traffic-light.png")
```

## 4. MFA analysis
MFA analysis done by `FactoMineR` library in R by F. Husson, J. Mazet algorythm implementation in sense of Esofier-Pages with supplementary individuals and groups of variables.

### 4.1. Factor proportions
There is explained variation, eigenvalue and cumulative percentage by factors and factor groups shown. 

```{r echo=FALSE}
tmp <- get_eigenvalue(fit) %>%
  round(., digits = 3)
colnames(tmp) <- c("EigenValue", "Variance", "Cumulative Variance")
knitr::kable(tmp, caption = "Table 4.1.1. Factor eigenvalues & explained variance proportions")
```

```{r echo=FALSE, fig.cap = "Fig. 4.1.1. Factor to dimensions explained variance"}
knitr::include_graphics("output/mfa/scree-plot.png")
```

Main value is Proportion of Variance show the proportion of explained variance in all data. For example, FAC1 proportion variance = 45 said that the 45% of data changes are explained by impact of FAC1 component.

### 4.2. Factor contributions
Here the factor group (PC) contribution for each one factor shown for main groups (1-3). This output can be used to any other assay.

```{r echo=FALSE}
fit.impact$coord %>%
  round(., digits = 3) %>%
  knitr::kable(., caption = "Table 4.2.1. Factor contribution coordinates in MFA")
```



```{r echo=FALSE, fig.cap="Fig. 4.2.1. Factor group contribution to dimensions 1-3"}
knitr::include_graphics("output/mfa/factor-contrib-group.png")
```

```{r echo=FALSE, fig.cap="Fig. 4.2.2. Factor individual contribution to dimensions 1-3"}
knitr::include_graphics("output/mfa/factor-contrib-var.png")
```

```{r echo=FALSE, fig.cap = "Fig. 4.2.3. Factor year variance contribution to dimensions 1-3"}
knitr::include_graphics("output/mfa/factor-contrib-year.png")
```

### 4.3. Weighted rotated matrix
In this section raw rotated standard deviation matrix to each one pc-dimension shown. Ofter this matrix used to explain factor changes year-by-year and can be used to explaine regime shifts.

```{r echo=FALSE}
fit$ind$coord %>%
  round(., digits = 3) %>%
  knitr::kable(., caption = "Table 4.3.1. Stdev DIM rotation matrix")
```

### 4.4. Main dimensions biplot
There is a 2-dimensional (biplot) scientific graphs explaining DIM1 changes over DIM2 (DIM1 - 1st dimension, DIM - 2nd dimension in xy plot). Vectors with same directions show the positive correlation between them. Opositve vectors show negative corelated factors.

```{r echo=FALSE, fig.cap="Fig. 4.4.1. Biplot DIM1 vs DIM2"}
knitr::include_graphics("output/mfa/biplot-dim1-vs-dim2.png")
```

```{r echo=FALSE, fig.cap = "Fig. 4.4.2. Individual biplot for years DIM1 vs DIM2"}
knitr::include_graphics("output/mfa/years-dim1-vs-dim2.png")
```

```{r echo=FALSE, testgl, webgl=TRUE, fig.cap="Fig 4.4.3. 3D plot DIM 1-3 year variance"}
library("rgl")
#knitr::knit_hooks$set(webgl = hook_webgl)
fit.coord <- fit.score$coord[,1:3]
fit.coord.offset <- fit.coord * 0.1
plot3d(x = fit.coord[,1], 
       y = fit.coord[,2], 
       z = fit.coord[,3],
       xlab = "Dim 1",
       ylab = "Dim 2",
       zlab = "Dim 3",
       size = 10, 
       col = "#0d0e69")
lines3d(x = fit.coord[,1], y = fit.coord[,2], z = fit.coord[,3], lty = "dashed")
text3d(x = fit.coord[,1] + fit.coord.offset[,1], 
       y = fit.coord[,2] + fit.coord.offset[,2], 
       z = fit.coord[,3] + fit.coord.offset[,3], 
       rownames(fit.coord), 
       col = "black",
       cex = 0.7,
)
rglwidget()
```

### 4.5. Years with trajectories
This plot allow to research trajectories of dimensions for single years to prove changes in time effect. 
```{r echo=FALSE, fig.cap = "Fig. 4.5.1. Individual biplot for years with trajectories DIM1-3"}
knitr::include_graphics("output/mfa/biplot-years-traj.png")
```

## 5. Regime shift detection
Regime shift detection approach based on PCA results. Instant changes in DIM stdevs show the group regime shifts. 

### 5.1. Dimension var changes

```{r echo=FALSE, fig.cap="Fig. 5.1.1. Dimension stdev center distance changes"}
knitr::include_graphics("output/mfa/preliminary-regime-shift.png")
```

### 5.2. Rodionov regime-shift
There is typical Rodeonov (2004) regime shift detection alorithm pefrormed. Mode: mean, p = 0.05, L = 10. Regime shift summary plot shown below. 

```{r echo=FALSE, fig.cap = "Fig. 5.2.1. Rodionov regime shift detection output"}
knitr::include_graphics("output/mfa/pca-regime-shift-radeonov.png")
```

### 5.3. Hierarchical clustering
In this section hierarchical cluster analysis for dissimilarities shown. Cluster analysis done by hclust() function with method “ward.D2”. Ward2 method use squared Euclidian distance minimization between component variances to determine hierarchical relations.

```{r echo=FALSE, fig.cap = "Fig. 5.3.1. Hierarchical cluster analysis for MFA variations"}
knitr::include_graphics("output/mfa/cluster-dendrogram.png")
```

## 6. Pearson correlation relationships
This section shows raw cross-correlation results for all data values. This analysis design is a typical for 1980-1990 years but not useful for present days. 

```{r echo=FALSE}
knitr::kable(cor.mat, digits = 2, caption = "Table 6.1. Pearson coefficient values (r) output")
```

In the figure below the summary short results shown. Only colored cells are significant and they correlation is true. 

```{r echo=FALSE, fig.cap="Fig. 6.1. Cross-correlation Pearson test on sign.level = 0.05"}
knitr::include_graphics("output/correlation-test-col.png")
```